function [bfinal1 pval mpy1 mpx mpu1 converg1] = fapr1(yy, ylag, x, opt, nsim, mmax, k, N)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                           %%
%%%% Factor Augmented Dynamic Panel Regression %%
%%%%  more than two variables                  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%-------------------------------------------------------------------
%%%
%%% yy(it) = a*ylag(it-1) + b*x(it) + c(i)*F(t) + e(it)
%%% yy - dependent variable (T x N)
%%% ylag - lagged dependent variable (T x N)
%%% x - explanatory variables (T x N(k-1)
%%% opt - factor selection criteria
%%% nsim - number of iteration for the modified PPC estimator
%%% mmax - maximum number of common factors 
%%% k - number of regressors including the lagged dependent variable, ylag
%%% N - number of cross-sectional units
%%% T - number of time periods 
%%%--------------------------------------------------------------------
%%%     
%%%     Date: June 3, 2013
%%%     Khusrav Gaibulloev, American University of Sharjah 
%%%     Todd Sandler, University of Texas at Dallas 
%%%     Donggyu Sul, Univeristy of Texas at Dallas
%%%---------------------------------------------------------------------


%%%%%%  Preparing data
xxx = [ylag x]; % a matrix of regressors

i=1; j=N; m=1;
Xcell=cell(1,k); % cell array of k empty matrices; will be filled below 

while (i<=((k-1)*N+1)) && (j<=k*N) && (m<=k)
    Xcell{m} = xxx(:,i:j); % each matrix in the cell represents a variable, TxN
    i=i+N;
    j=j+N;
    m=m+1;
end


%%%%%%  Withing Group estimation: residuals are used for estimating the
%%%%%%  number of common factors next
for r=1:k
    xx(:,r)=Xcell{r}(:);
end
beta1 = (xx'*xx)\(xx'*yy(:));

xx=[];
for i=1:N
    for r=1:k
        xxtemp(:,r)=Xcell{r}(:,i);
    end
    xx=[xx xxtemp];
end

clear xxtemp;
beta1=kron(eye(N),beta1);
res1 = yy - xx*beta1;


%%%%%%  1. Estimating the number of common factors
mpy1 = ppbn(yy,mmax,opt,1); %TxN, max common factors, criteria; prewhitening=1;
mpx = ppbn(xxx, mmax,opt,1);
mpu1 = ppbn(res1,mpy1,opt,0);


%%%%%%  2.  The PPC estimation  

%%%%%%  2a. Estimating common factors from regressand and regressors
[fy1,ly1,my1,vnt] = pce(yy,mpy1);
[fx,lx,mx,vnt] = pce(xxx,mpx);

fac1 = [ones(size(fx,1),1) fy1 fx]; 
fqq1 = fac1*inv(fac1'*fac1)*fac1';

%%%%%%  2b. Obtaining the PPC Estimator
yf1 = yy;
xf1 = cell(1,k);

for i = 1:N 
    yf1(:,i) = yy(:,i) - fqq1*yy(:,i); 
        for j=1:k
            xf1{j}(:,i)=Xcell{j}(:,i)-fqq1*Xcell{j}(:,i);
        end
end

xx1=[];
for r=1:k
    xx1(:,r)=xf1{r}(:);
end

betaf1 = (xx1'*xx1)\(xx1'*yf1(:));


%%%%%%  3. The modified PPC estimation

%%%%%%  3a. Iterate until convergence: IPC using the PPC estimates as initial values 

yl1=yy;
xl1=cell(1,k);
xx1=[];
xx=[];

for i=1:N
    for r=1:k
        xxtemp(:,r)=Xcell{r}(:,i);
    end
    xx=[xx xxtemp];
end

clear xxtemp;
betaf = zeros(k,nsim);
betaf(:,1) = betaf1;
%betaf(:,1)=kron(eye(N),betaf);
s = 2;

while s <= nsim
    betaftemp = kron(eye(N),betaf(:,s-1));
    %res1s = yy - ylag*betaf(1,s-1) - tlag*betaf(2,s-1) - ws*betaf(3,s-1) - wn*betaf(4,s-1);
    res1s = yy - xx*betaftemp;
    [fu1,lcx1,mx1,vnt] = pce(res1s,mpu1);
    fau1 = [ones(size(fy1,1),1) fu1];           
    fxx  = fau1*inv(fau1'*fau1)*fau1';
        for i = 1:N
            yl1(:,i) = yy(:,i) - fxx*yy(:,i);
                for j=1:k
                    xl1{j}(:,i)=Xcell{j}(:,i)-fxx*Xcell{j}(:,i);
                end
        end
                        for r=1:k
                            xx1(:,r)=xl1{r}(:);
                        end
    betaf(:,s) = (xx1'*xx1)\(xx1'*yl1(:));
    if abs(betaf(:,s) - betaf(:,s-1))<0.0001; converg1 = 1; bfinal1 = betaf(:,s); break; 
    else converg1 = 0; bfinal1 = 'NA'; s = s+1; end;
end

%%%%%   3b. Standard errors and t-stats
bb = [];
xxx1=[];

for i = 1:N
    for r=1:k
        xxx1(:,r)=xl1{r}(:,i);
    end
        bb1 = (xxx1'*xxx1)\(xxx1'*yl1(:,i));
        bb = [bb bb1];
end

bb = bb';
bb = bb - repmat(mean(bb),N,1);
bb = bb'*bb/N;

for r=1:k
    xx1(:,r)=xl1{r}(:);
end
[t,n] = size(xl1{1});

rst1 = (xx1'*xx1).*bb.*(xx1'*xx1);
rst1 = rst1/(n*n*(n-1)*t*t);
rst11 = (xx1'*xx1)/(n*t);
rst11 = (rst11.*rst11)\rst1;
rst11 = rst11*n;
tt1 = (bfinal1)./sqrt(diag(rst11));

pval = 2*(1-normcdf(abs(tt1),0,1)); %compute p-value from the standard normal 


